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Abstract 

We investigate the Riemann problem for the shallow water equations with variable and (possibly) 
discontinuous topography and provide a complete description of the properties of its solutions: 
existence; uniqueness in the non-resonant regime; multiple solutions in the resonant regime. This 
analysis leads us to a numerical algorithm that provides one with a Riemann solver. Next, we 
introduce a Godunov-type scheme based on this Riemann solver, which is well-balanced and 
of quasi-conservative form. Finally, we present numerical experiments which demonstrate the 
convergence of the proposed scheme even in the resonance regime, except in the limiting situation 
when Riemann data precisely belong to the resonance hypersurface. 

Keywords: Shallow water model, hyperbolic conservation law, discontinuous topography, 
resonant regime, Riemann solver, Godunov-type scheme. 



1. Introduction 



In this paper we design a Godunov-type scheme for the numerical approximation of weak 
solutions to the initial-value problem associated with the shallow water equations with variable 
topography, i.e. 

dth + d.,{hu)^Q, 



dt{hu) + d^(h{u'^ + S^)^ +ghd.^a = 0, 



where the height of the water from the bottom to the surface, denoted by h, and the fluid velocity 
u are the main unknowns. Here, g is the so-called gravity constant, and a — a{x) (with x £ R) is 
the height of the bottom from a given level. 



In [27^, LeFloch pointed out that by supplementing balance laws like (1.1 ) with the additional 
equation 

dta = 0, (1.2) 



the set of equations ( 1.1 1-( 1.2 ) can be regarded as a nonlinear hyperbolic system in nonconservative 
form and tackled with the theory introduced by Dal Maso, LeFloch, and Murat [H] and developed 



by LeFloch and collaborators [23 [Ml [HI HO] ■ As is well-known, the system ( 1.1 )-( 1.2 ) is hyperbolic 
but not strictly hyperbolic since characteristic speeds may coincide on certain hypersurfaces. Non- 
strictly hyperbolic systems have been extensively studied in the literature. See [371 [Ml [Ml I3Z] for 
the model of fluid flows in a nozzle with variable cross-section, [TSl [31 S] for other models. 

On the other hand, the discretization of source terms in nonlinear hyperbolic balance laws 
like the shallow water equations was pioneered by Greenberg and Leroux [T^. We built here on 
this paper as well as the follow-up work [Tl[9l[8l[Ill[16lll8l|22l[23. Since the system (O) is 



also related to the class of two-phase models, we are also motivated by the existing research work 
on the discretization of two-phase flow models |5J[31[T3]. In particular, recall that well-balanced 
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schemes for shallow water equations were constructed first in [HI [171 1221 [23]. In addition, the 
discretization of nonconservative hyperbolic systems and of systems with source terms attracted a 
lot of attention in recent years. We refer to ^^jfl [TB] [IB] for a single conservation law with source 
term and to [531 [221 1251 121] for fluid flows in a nozzle with variable cross-section. Well-balanced 
schemes for multi-phase flows and other models were studied in [2 [3 [3H1 [3H] . 

The Godunov scheme is based on an exact or approximate Riemann solver and, consequently, it 
is necessary that sufficient information be available on the existence and properties of all solutions 



to (1.1|-(1.3|. Recall that the Riemann problem is a Cauchy problem with piecewise constant 



initial data of the form 

(h,u,a)(x,0) = < ;, \ „ (1.3) 

One main task in the present paper will be to revisit the construction of solutions to the Riemann 



problem for (1.1)-(1.3) in order to arrive at a definite algorithm for their numerical computation. 
We will show that a unique solution exists within a large domain of initial data, and will precisely 
identify the domains where multiple solutions are available, by providing necessary and sufficient 
conditions for the existence of multiple (up to three) solutions. 

Recall that, in LeFloch and Thanh [53], a first investigation of the Riemann problem for the 
shallow water equations was performed. The present paper provides a very significant improvement 
in that a complete description of a Riemann solver is now obtained for the first time. We are able 
to distinguish between cases of existence, uniqueness, and multiplicity of solutions. 

We refer the reader to [U [HI [35] for partial or alternative approaches to the Riemann problem. 



On the other hand, Chinnayya, LeRoux, and Seguin [TTl introduced a Godunov method for (1.1 1 
based on a Riemann solver determined by "continuation" : they start their construction by assum- 
ing that the bottom topography is flat and then extend it to a non-flat topography. Their method 
allows them to construct solutions within the regime where the nonlinear characteristic fields are 
separated by the linearly degenerate one, that is, the case where one wave speed is negative and the 
other positive. In this regime, the Riemann solution with non-flat bottom can be obtained from 
the the wave curves associated with the fast wave family, only, and these two waves are separated 
by a stationary wave. In particular, the total number of waves in Riemann solutions is exactly 
the number of characteristic families. On the other hand, when more general Riemann data are 
considered and lie in the other two strictly hyperbolic regions, say of "cross-region" type, we find 
it hard to apply this method. There is always a major difference between the "fiat-bottom" and 
"non-flat-bottom" cases, since in one case the system is strictly hyperbolic, the other case the 
system is not strictly hyperbolic. Indeed, in the non-flat-bottom case, new wave curves arise that 
replace more standard wave curves. The total number of waves in a solution can possibly be larger 
than the number of characteristic flelds as waves associated with a given family can be repeated. 
The appearance of such new wave curves cannot be obtained by a continuation argument. 

Our objective in the present work is to present a numerical algorithm that provides an explicit 
construction of a Riemann solver for the shallow water equations in the resonant or non-resonant 
regimes. This solver then provides us with solutions to local Riemann problems which we can use 
to design a Godunov-type scheme. We also provide extensive numerical experiments and, within 
strictly hyperbolic regions of the phase space, our tests demonstrate that the proposed scheme 
converge to the expected solution. In resonance regions, we also observe convergence, except 
when the Godunov scheme takes some values on the resonance hyper surf aces. The Godunov 



scheme proposed in the present paper for (1.1 1 (1.21 turns out to be well-balanced and captures 
exactly stationary waves. 



This paper is organized as follows. In Section 2, we recall basic facts about the system (1.1) 



( 1.2 ) and provide the computing algorithm for stationary contact waves. In Section 3, we revisit the 
construction of solutions to the Riemann problem, and present a new approach based on a "gluing" 
technique involving different solution structures. With this technique, we establish the existence 
of solutions for arbitrary large initial data. In Section 4 we then present our computing strategy 
leading to a Riemann solver and, finally, are in a position to design a corresponding Godunov 
scheme. Section 5 is devoted to numerical experiments when data belong to strictly hyperbolic 
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domains; in particular, we estimate the numerical errors and observe convergence when the mesh 
size tends to zero. Finally, Section 6 is devoted to numerical tests here data belong to resonance 
regions, and precisely identify regimes of convergence or non-convergence. 



2. Shallow water equations 
2.1. Wave curves 

Introducing the dependent variable {h, u, a) — {h, u, a)(x, t), the Jacobian matrix of the system 
(1.1)-(1.2| admits three real eigenvalues, i.e. 

Xi{U):=u-y^<X2{U):=u+^, X3{U):=0, (2.1) 

together with the corresponding eigenvectors: 

gh 

n{U) := \ -Vgh \ , r2(t/) := I VP I , r,iU):=\ -gu I, (2.2) 





,2 



so that the system (1.1)-(1.2) is hyperbolic, but not strictly hyperbolic. More precisely, the first 
and the third characteristic speeds coincide, 

(Ai(C/),ri(C/)) = (Z3(C/),r3(C/)), 

on the hypersurface 

C+~[{Ku,a)\ u-yP}. (2.3) 
The second and the third characteristic fields coincide, 

(A2(C/),r2(C/)) = (A3(t/),r3(C/)), 

on the hypersurface 



C_ := 



[{h,u,a)\ u = -y/^y (2.4) 



In the following, we introduce C := C_|_ U C_, and we refer to the sets C± as the resonance 
hypersurfaces, which separate the phase space in the {h, u, a)-variable into three sub-domains 



Gi := Uh,u,a) e R+ x R x R+\ Xi{U) > XsiU)^ 



G2 



ih,u,a) eR+xRxR+\ A2(C/) > A3 (t/) > Ai ([/)}, (2.5) 
G3 := {{h,u,a) eR+xRxR+\ XsiU) > A2([/)}, 

in which the system is strictly hyperbolic. It is convenient to also set 

Gt {{h, u, a) e G2, M > 0}, {{h, w, a) e G2, u < 0}. 

Observe that for u > 0, the set Gi is referred to as the domain of supercritical flows, where the 
Froude number 

„ u 
Fr := ^= 

is larger than 1, and the set Gj as the domain of subcritical flows, where Fr < 1. This concept 
may similarly be extended to the case u < 0. See Figure [ij 



One easily checked that the first and second characteristic fields (Ai,ri), (A2,r2) are genuinely 
nonlinear, while the third characteristic field (A3,r3) is linearly degenerate. 
As discussed in [33], across a discontinuity there are two possibilities: 
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Figure 1: Phase domain in the (h,M)-plane. 



(i) either the bottom height a remains constant, 

(ii) or the discontinuity is stationary (i.e. propagates with zero speed). 



In the first the case (i), the system 1.2 1 reduces to the standard shallow water equations 



with flat bottom. We can determine the i-shock curve Si{Uo) starting from a left-hand state Uq 
and consisting of all right-hand states U that can be connected to Uq by a Lax shock associated 
with the first characteristic field (i = 1,2): 



5,(C/o): ^,{U;Uo):^u~uo±^ih-ho)^lQ~+^ = 0. (2.6) 

where U = {h,u), h > Hq for i = 1 and h < Hq for i = 2. Wc also define the backward j-shock 
curve SP{Uo) starting from a right-hand state Ua and consisting of all left-hand states U that can 
be connected to Uq by a Lax shock associated with the first characteristic field {i — 1,2): 

5f(C/o): MU; Uo) :=u-uo± ^-{h - h^)^ (1 + 1) = 0, (2.7) 

where U = {h, u), h < ho for i — 1 and h > ho for i ~ 2. 

It is interesting that the shock speed of the nonlinear characteristic fields may coincide with 
the speed of stationary contact waves. The following lemma is easily checked. 

Lemma 2.1. Consider the projection on the (h,u)-plan. To every U — {h,u) G Gi there exists 
exactly one point U"^ £ Si{U) H such that the 1-shock speed Xi{U,U'^) = 0. The state 
U"^ = {h'^^w^) is defined by 

2 ' " =X#- 

Moreover, for any V G Si{U), the shock speed Xi{U, V) > if and only if V is located above t/^ 
on Si{U). 

It is also well-known that the bottom height a remains constant through rarefaction fans. The 
forward rarefaction curve TZi{Uo) starting from a given left-hand state Uq and consisting of all the 
right-hand states U that can be connected to Uq by a rarefaction wave associate with the first 
characteristic field as 

7^,(C/o): ^^{U;Uo)^u-uo±2^{Vh-^/h^)^0, ^ = l,2, (2.8) 
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where U = {h,u),h < ho for i = I and h > ho for i = 2. Given a right-hand state Uq, the 
backward i-rarefaction curve TZf{Uo) consisting of all left-hand states U that can be connected to 
Uq by a rarefaction wave associated with the first characteristic field reads (i = 1, 2) 



7^f(^7o): MU;Uo)^u~UQ±2^{Vh~^o)^0, (2.9) 

where U = {h,u), h > ho for i — 1 and h < ho for i = 2. 

Finally, we define the forward and backward wave curves in the (/i,M)-plane (z = 1,2): 



W.(C/o) := S,{Uo) U 7^,([/o) - {U \ ^,{U; Uo) = 0}, 

Wf (C/o) := 5f (C/o) U 7^f (C/o) = {U \ ^U; Uo) = 0} ^^'^^^ 



It is checked in [33^ that the wave curves Wi(C/o) and yVi{Uo) parameterized as h i-^ u — u{h), h > 
0, are strictly convex and strictly decreasing functions. The wave curve yV2(t^o) and yV^(?7o) being 
parameterized as /i i— > u = u{h), h > 0, are strictly concave and strictly decreasing functions. 
In the case (ii), the discontinuity satisfies the jump relations 



(2.11) 



[hu] = 0, 

[y^+g{h + a)]=0, 

which determine the stationary- wave curve (parameterized with h): 

mUo): u^uih) = ^,^^ ^^^^^^ 
a = a{h) = flo -f- "°2g" — I" '^0 ^ h. 

The projection of the wave curve W^^Uo) in the (/i, u)-plane can be parameterized as h t-^ u = 
u{h),h > 0, which is a strictly convex and strictly decreasing function for uq > and strictly 
concave and strictly increasing function for uq < 0. 

The above arguments show that the a-component of Riemann solutions may change only across 
a stationary wave. This property will be important later when designing the discretization of the 
source terms. 

2.2. Properties of stationary contacts 

Given a state Uq = {ho,uo,ao) and another bottom level a ^ oq, we let U — {h,u,a) be the 
corresponding right-hand state of the stationary contact issuing from the given left-hand state 
Uo- We now determine h,u in terms of Uo^a, as follows. Substituting u = houo/h from the first 



equation of (2.121 to the second equation of (2.12), we obtain 

2\ 




/ hpUQ \ ' 



ho — h = a. 



Multiplying both sides of the last equation by 2gh^, and then re-arranging terms, we find that 
/i > is a root of the nonlinear equation 

ip{h)^Lp{Uo,a\h) 

:=2gh^ + {2g{a-ao-ho)-ul)h^ + hlul^Q. 



We easily check 



so that 



(^(0) - hlul > 0, 

ip'{h) = Qgh"^ + 2{2g{a - - ho) - ul)h, 
ip"{h) = Ugh + 2{2g{a - ao - ho) - u^), 

ip'{h) = iff h = 

ul + 2g{ao + ho - a) (2-14) 



or h = h^ = h^,{Uo,a) := 
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Figure 2: Graph of the function (p = fih), h > defined by l |2.13[ l with g = 9.8, ao = l,hQ = 1,uq = 1 and a = 1.2. 
The function ifi admits two zeros in the interval (0, 1). 



If h^,{Uo, a) < 0, or a > flp + /lo + then ip'{h) > for h > 0. Since tp{0) — h^ul > 0, there is 



no root for (2.13) if (2.14) holds. Otherwise, if 



a<ao + ho + —, 

then (fi' > ior h > ft,* and (p'{h) < for < /i < /i*. In this case, ip admits a zero h > 0, and in 
this case it has two zeros, iff 



>f{K) = -ghl + hlul < 0, 



where /i* is defined by (2.14). It is easy to check that (2.15) holds if and only if 



a < amax(C/o) := ao + /lo + 1^ - 5^173 (/louo)^/^ 



ao 



2g 



(2.15) 



(2.16) 



Observe that (2.16) implies aniax(C^o) ^ and the equality holds only if {hQ,UQ) belongs to the 



surfaces C±. Whenever (2.16) is fulfilled, the function ip in (2.131 admits two roots denoted by 
hi{a) < h2ia) satisfying hi{a) < /i* < /12(a). Moreover, if the inequality in ( |2.16 ) is strict, i.e., 
a < amax(t^o), then these two roots are distinct: hi{a) < /i* < /12(a). See Figure [2j 



Thus, we arrive at the following lemma. 
Lemma 2.2. Given Uq ~ {Hq, uq, gq) and a bottom level a =/= aQ. The following conclusions holds, 
(i) amax(C/o) > ao, amax(C/o) = ao if and only if {Ho^uq) G C. 



(ii) The nonlinear equation (2.13 ) admits a root if and only if the condition (2.16) holds, and in 
this case it has two roots /ii(a) < /i* < /12(a). Moreover, whenever the inequality in (2.16) 
is strict, i.e. a < aniax(t^); these two roots are distinct. 
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(in) According to the part (ii) , whenever (2.16) is fulfilled, there are two states Ui{a) = [hi{a),Ui{a),a), 
where Ui{a) = hQUQ/hi{a),i = 1,2 to which a stationary contact from Uq is possible. More- 
over, the locations of these states can be determined as follows 



Ui{a) € Gi 
Ui{a) e G3 
f/2(a) € G2. 



Uq > 0, 
Uo < 0, 



Proof. The parts (i) and (ii) can be easily deduced from the above argument. To prove (iii), it is 
sufficient to show that along the projection of W3{Uo) on the {h,u)-p\anc, the point J7min(C^o) — 
(/imin(C^o),'"min(t^o) ^oWq/ Vin (C^o) ) , where /imin(t^o) is defined by ( |2.15[ ), belongs to C+ if 
Mo > and belongs to C_ if uq < 0, and that Ui{a) € yV3(C/o),i = 1,2, such that C/2(a) £ G2 
and t7i(a) is located on the other side of [/2(a) with respect to C. Indeed, let us define a function 
taking values along the stationary curve 'W3{Uo): 



a(h) -.^ ulhf - gh = ^ - gh. 

Clearly, a point U — {h,u,a) belongs to Gi U G3 if and only if a{h) > and U belongs to G2 
if and only if (T{h) < 0. Since a-{h^in{Uo j) — 0, the point UminiUo) belongs to C. Obviously, 

(Uq) E C_ if Uq < 0. Thus, it remains to check that 



Umin{Uo) e C+ if Uo > 0, and C/, 



cr(/ii(o)) > 0, a{h2{a))<0. 



(2.17) 



Since 



a(/imi„(C/ - 0)) - 0, <j'ih) = 



-5<0, 



we can see that (2.17) holds if 



hi{a) < hmin{Uo) < h2{a). 



On the other hand, we have 



(p[h) > 0, iih < hi{a) or h > /12(a), 
(f{h) < 0, if hi{a) < h < /12(a). 



(2.18) 



(2.19) 



And we have 



^(hminiUa)) = 3(/ioUo)^ + (2.9(a - ao - ho) - ul) 



72/3 



It is a straightforward calculation to show that the condition 



a < a„ 



is equivalent to 



</'(/lmin(C/o)) < 0. 



This together with (|2.19|) establish (|2.18|). Lemma 2.2 is completely proved. □ 

The Riemann problem may 



From Lemma |2.2[ we can construct two-parameter wave sets, 
therefore admit up to a one-parameter family of solutions. To select a unique solution, wc impose 
an admissibility condition for stationary contacts, referred to as the Monotonicity Criterion and 
defined as follows: 

(MC) Along any stationary curve W3(C/o), the bottom level a is monotone as a function of h. The 
total variation of the bottom level component of any Riemann solution must not exceed 
\o-L ~ o-r\: where a^, a^j are left-hand and right-hand bottom levels. 



A similar criterion was used [201 EI 



and \m 
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Lemma 2.3. The Monotonicity Criterion implies that any stationary shock does not cross the 
boundary of strict hyperholicity, in other words: 

(i) If Uq S Gi U G3, then only the stationary contact based on the value hi{a) is allowed, and 
one sets h{a) — hi{a). 

(a) If Uq £ G2, then only the stationary contact using h2{a) is allowed, and one sets h(a) — 
/12(a). 

Thus, h{a) is the admissible /i- value of a right-hand state U — {h — h{a), u, a) of the stationary 
wave from a given left-hand state Uq = (/io,uo,ao). 

Proof. Recall that the Rankine-Hugoniot relations associate the linearly degenerate field (2.11 1 
implies that the component a can be expressed as a function of h: 

a = a{h) = ao H n + ho, 

where 

u — u(h) = . 

h 

Thus, taking the derivative of a with respect to h, we have 

a'{h) ^^^-l^u^f^-l 

gh gh 

which has the same sign as — gh. Thus, a = a{h) is increasing with respect to h in the domains 
Gi,G3 and is decreasing in the domain G2. Thus, in order that a = a{h) is monotone as a 
function of h, the point {h,u,a) must stay in the closure of the domain containing (hQ,UQ,ao). 
The conclusions of (i) and (ii) then follow. □ 
We now explain how to compute the roots of the equation (2.13). The above argument shows 
that whenever (2.16) is satisfied, the equation (2.13) admits two roots /ii(a), /12(a) satisfying 

h,(a) < /._ = C^y^' < K - + + < h.ia) (2.20) 



and the inequalities are all strict whenever the inequality in (2.16) is strict. Since < hi{a) < 
hmin < ^* and 

^(0) > 0, 

<fi{K) < 0, (^(/imin) < 0, 



the root /ii(a) of (2.13) can be computed, for instance using the regula falsi method with the 
starting interval [0,/imin], or [0,/i*]. And since /12(a) > /i* and (p'{h) > 0,ip"{h) > 0,h> /i*, the 
root /12(a) can be computed using Newton's method with any starting point larger than /i*. We 
summarize this in the following lemma. 



Proposition 2.4 (Water height of stationary contacts). The root hi{a) of {2.13) can be computed 

(h^u^ \ ^^'^ 
"J'" j , or [0, h^,], 

where h^ — "^o+^siao+ho a) ^ ^^j^^ ^j^g /12(a) can be computed using Newton's method with any 
starting point larger than /i* . 

To conclude this section, we point out that certain physical applications may actually require a 
different jump relation for the nonconservative product — especially allowing for energy dissipation. 
This issue will not be addressed further in the present paper, however. 



8 



3. The Riemann problem revisited 

From the general theory of nonconservative systems of balance laws, it is known that if Riemann 
data belongs to a sufficiently small ball in a strictly hyperbolic region, then the Riemann problem 
admits a unique solution. It is worth to note that this result no longer holds if any of these 
assumptions fails, for instance due to resonance. 

Our goal in this section to to provide all possible explicit constructions for Riemann solutions, 
investigating when data are around the strictly hyperbolic boundary C±. There are several im- 
provements in the construction of Riemann solutions in this paper over the ones in our previous 
work '333- First, we can determine larger domains of existence by combining constructions in [3 3) 
together. Second, the domains where there is a unique solution or there are several solutions are 
precisely determined. Under the transformation x i— —x,u i— > —u, a left-hand state U — {h,u,a) 
in G2 or G3 will be transferred to the right-hand state V = {h, —u,a) in G2 or Gi, respectively. 
Thus, the construction for Riemann data around C_ can be obtained from the one for Riemann 
data around C+ . We thus construct only the case where Riemann data are in Gi U C+ U G2 and 
we separate into two regimes on which a corresponding construction based on the left-hand state 
Ul is given: 

• Regime (A): Ul e Gi UC+; 

• Regime (B): Ul G G2; 

For each construction, depending on the location of the right-hand states Un and the sign a/? — 
there will be different types of solutions or the results on the existence and uniqueness. 
As in to solve ( |l.l[ )-( |L3l ) we project all the wave curves on the (/i, M)-plane. 



Notations 

(i) U^ denotes the state resulted from a stationary contact wave from U] 



(ii) U'^ is the state defined in Lemma 2.1 so that Ai(C/, [/*) — 



(hi) Wk{Ui,Uj) {Sk{Ui,Uj), Rk{Ui,Uj)) denotes the fcth-wave (fcth-shock, fcth-rarefaction wave, 
respectively) connecting the left-hand state Ui to the right-hand state Uj, k = 1,2,3; 

(iv) WmiUi,Uj) © Wn{Uj,Uk) indicates that there is an mth-wave from the left-hand state Ui 
to the right-hand state Uj, followed by an nth- wave from the left-hand state Uj to the 
right-hand state U^, m,n £ {1,2,3}. 

3.1. Regime (A). Eigenvalues at Ul with coinciding signs 

Let Gi denote the closure of Gi. We assume that Ul € Gi, or equivalently Xi{UL) > 0,i = 
1,2,3. 

Construction (Al). In this case (the projection on (/i, u)-plane of) Ur is located in a "higher" 
region containing Ul in the (ft,, M)-plane. 

If fli > Or (or aL < aji < a-^aai^iU l)) , the solution begins with a stationary contact upward 
(downward, respectively) along yV^{UL) from Ul to the state J7£ € yV3(C/L)nGi, shifting the level 
ttL directly to the level or. Let 

{Um = {hM,UM,aR)} - Wi{Ul) n Wi{UR). 



Providing that Ai(J7£, Um) > 0, or equivalently, as seen from Lemma 2.1 ftj\/ < h°j^ , the solution 
can continue by a 1-wave from Ul to Um, followed by a 2-wave from Um to Ur. Thus, the solution 
is 

W^{Ul, Ul) ® VFi([/£, Um) ® W2{Um, Ur). (3.1) 
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Figure 3: Construction (Al), > a.R'- a solution of the form | |3.1| l. 



See Figure [3] This construction can be extended if yV:f(?7fl) lies entirely above Wi(C/£). In 
this case let / and J be the intersection points of WiiUl) and >V^(/7fl) with the axis {h ~ 0}, 
respectively: 

{/} = Wi((7£)n{/i = o}, {j}^wi{UR)n{h^o}, (3.2) 

then the solution can be seen as a dry part Wo{I, J) between / and J. Thus, the solution in this 
case is 

W3iUL,Ul) ® RiiUlJ) © Wo{I, J) ® i?2(J, Ur). 



Remark 1. As seen by Lemma [2^2] if < a/j, the condition 

is necessary for the stationary contact W3(J7i,, C/£). Therefore, if this condition fails, there is no 
solut ion even it Ul = Ur. The necessary and sufficient conditions for the existence of the solution 
(3.1) is that C/^'^ is located below or on the curve {Ur). This domain clearly covers a large 



crossing-strictly- hyperbolic-boundary neighborhood of Ul. 

Construction (A2). In this construction, we will see an interesting phenomenon when wave speeds 
associate with different characteristic fields coincide. Roughly speaking, this case concerns with 
the fact that Ur moves limitedly downward from the case Gi. Instead of using "complete" 
stationary contact from Ul to C/£ as in the first possibility, the solution now begins with a "half- 
way" stationary contact Wz{Ul,Ui) from Ul — {h,u,aL) to some state Ui = f/£(a) = {h,u,a) £ 
W3(C/l), where a between ul and or. The solution then continues by a 1-shock wave with zero 
speed from Ui to U2 = Uf G Wi(J7i) n G2. Observe that U2 still belongs to WziUL), since 



/12W2 = hiui — Hlul, as indicated by Lemma 2.1 The solution continues by a stationary contact 



from U2 to a state Um{<i) € W3(C/l). The set of these points UM{o),a £ [aL,aR\ forms a curve 
pattern denoted by £. Whenever 

Wf(i7fl,)n£^0 

there is a solution containing three discontinuities having the same zero speed of the form 

WsiUL, Ul) © SiiUi,U2) © W3iU2,UM) © W2{Um, Ur). (3.3) 

See Figure |4] 
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Figure 4: Construction (A2), ax, > an: a solution of the form | |3.3| l. 

Remark 2. The necessary and sufficient conditions for the existence of the solution ( |3.3[ ) is that 
U'^ is located above or on the curve Wf (C/fl), and U*° is located below or on the curve 'W2 {Uji). 
This domain covers a region in and Gi which is "far away" from U^. 

It is interesting that at the limit a = aj^ at the first jump, we get the first possibility. If 
a — ai,, then the solution simply begins with a 1-shock wave with zero speed followed by a 
stationary contact shifting a from to ur. This limit case can be connected to the following 
possibility. 

Construction (A3). The solution begins with a strong 1-shock wave from Ul to any state U G 
Wi(C/l) n G2 such that Xi{Ul, U) < 0. This shock wave is followed by a stationary contact to a 
state U° shifting a from to an. The set of these states U° form a curve denoted by (?7l)- 
That is 

((7l) {[/° : 3W3{U,U°) shifting az. to a^,, 

U^{h,u,aL)eWi{UL),XiiU)<0}. > 

Whenever 

0^Wf(C/fl)nWi"(f/L) = {C/^JcG2, and MUIi^Ur) > 0, (3.5) 
there will be a Riemann solution defined by 

Si{Ul, Um) © W3{Um, UIj) ® W2{Uli, Ur). (3.6) 



See Figu re [Sj In the limit case of (3.3) where Ui = Ul, the solution (3.3) coincides with the 



solution (3.6) 



Let K denote the lower limit state on Wi{Ul) that the solution (3.6 1 makes sense, and let 
K" £ G2 denote the right-hand state resulted from a stationary contact from K shifting ol to an. 
Thus, we have 

{X} = Wi(C/z,)nC_, ifai>afl, 

KeWi{UL) such that ai„ax(-?f) = a_R, if < a_R. 



Remark 3. The solution (3.6) makes sense if Uf° is above or on the curve yV-^iUn), and K° lies 



below or on the curve (CTfl) and A3 (if", Ur) > 0. 
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Figure 5: Construction (A3), aj^ > ajj: a solution of the form ||3.6| 



The union of the wave patterns VVi(C/i) U £ U W°(C/l) form a continuous curve. The Riemann 
problem thus admits a solution whenever WHUr) intersects WiiUi) UCU Wi{Ul) or WHUr) 
intersects with {h = 0} at a point above the point /. We can see that this happens for a large 
domain of Un containing Ul ■ See Figures [3j [4j and [5j 

If the wave pattern £ lies entirely on one side with respect to the curve WfiUB), then yV;f (?7ij) 



intersects either >Vi(f/i) or Wi(Ul) at most one point. Therefore, then (3.1) or (3.6 1 is the unique 



U 



■)lution. Besides, if yV^(C/fl) intersects the wave pattern £, and if hf^° > h°j* , then the point 

on the curve W3(C/l). Thus, the curve W^(J7ij) does not 



^ is located below the point Uj^ 



meet Wi{Ui) nor >Vi(t7i,), except possibly at the endpoints Uf° € C and t/^'^ G C. In this case, 
(3.3) is the sole solution. In summary, the Riemann problem for ( |l.l[ )-( [L2l ) always has at most 



one solution whenever h 



In the case where hf^° < hj^ 



there can be three solutions, as yV^iUn) can meet all the three 
curve patterns yVi{Ui),£ and (C/l), or 



< 



h°* 

"■L ' 



(3.8) 



where the function ^2iU; Ur) is defined by (2.10). See Figure|6] 
The above argument leads us to the following theorem. 

Theorem 3.1 (Riemann problem for the shallow water equations). Given a left-hand state Ul ^ 
Gi . Depending on the location of the right-hand state Ur we have the following conclusions. 



(a) Existence. The Riemann problem (1.1 )- ^T7§^ admits a solution if Q° defined in Construc- 
tion (A3) lies below or on the curve WiiUii), and that ifW^iUn) intersects with W^{Ul) 
at some point U^j G G2 , then A2 (?/£/, Ur) > 0. 



(b) Regime of uniqueness. The Riemann problem (1.1)- (1.3) has at most one solution if 



- either h*° > h°* ; 

- or < , and the states Uf(° and U'^ are located on the same side with respect to 
the curve W|^(C/k). 



12 



Figure 6: Non-uniqueness: three admissible Riemann solutions of the form l |3.1| l, ( |3.3| l, and | |3.6[ |. 
(c) Multiple solutions. If hf^° < h"^ , and if the state Uf" lies above the curve (Ur) 



while the state lies below the curve 'W2{Ur), then the Riemann problem {1.1)-(1.3) 
has three solutions. 

Example. We provide some numerical experiments to illustrate two situations: > h"j^ , and 
hf° < h°* corresponding to the two cases ql > clr (see Tables A1-A3) and < (see Tables 
A4, A5). We take at random the state Ul and a_R. 

(a) > a^: all experiments show that > h°j^ . In Table Al, G Gi. 



Table Al 


States 


Ul 


ur 


uf 


Water Height h 


0.5 


1.1930011 


1.1171275 


Velocity u 


4 


1.6764444 


1.790306 


Bottom Level a 


1 


0.9 


0.9 



In Table A2, Ul e 



Table A2 


States 


Ul 


ur 


uf 


Water Height h 
Velocity u 
Bottom Level a 


1 

3.1304952 
1 


1.3075478 
2.3941726 
0.9 


1.2558035 
2.4928225 
0.9 



In Table A3, [/l G Gi is far away from C+. 



Table A3 



States 


Ul 


ur 


Ul* 


Water Height h 


0.01 


0.54763636 


0.44902891 


Velocity u 


10 


0.18260292 


0.22270281 


Bottom Level a 


1 


0.9 


0.9 
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(b) gl < aji'. all experiments show that hf° < . In Table A4, Ul E Gi. 



Table A4 


States 


Ul 


ur 


Uf 


Water Height h 


0.5 


0.86127059 


0.96534766 


Velocity u 


4 


2.3221506 


2.0717925 


Bottom Level a 


0.9 


1 


1 



In Table A5, [/l G Gi is far away from C+. 



Table A5 



States 


Ul 


ur 


uf 


Water Height h 


0.01 


1.2748668 


1.3718425 


Velocity u 


10 


0.78439566 


0.72894668 


Bottom Level a 


0.9 


1 


1 



Remark 4. We conjecture that if > o-r, then hf > hf , and if < an, then hf < hf . If 
this conjecture holds, then Theorem |3 . 1 1 implies that when > fli?, the Riemann problem has at 
most one solution for Ul S Gi. 

3.2. Regime (B). Eigenvalues at Ul with opposite signs 

In this subsection we consider the case where the left-hand state Ul moves downward from the 
Regime (A): Ul € G2, or Ai(C/l) < = MUl) < MUl)- 

Construction (Bl). For Ur in a "higher" position, there can be two types of solutions depending 
on whether > a^. 

li ol > oiR. a solution can be constructed as follows. The solution begins from Ul with a 
1-rarefaction wave until it reaches C^. at a state C/i G Wi(?7l) nC_|_. A straightforward calculation 
gives 

U^-[{^g + -i^hL) ,-^UL + -^^9hL,aL). 

This rarefaction wave can be followed by a stationary jump W3(?7i, U2) into Gi. This stationary 
wave is possible since ol > an. Let {C/3} = 'Wi{U2) H 'W2{Ur). The solution is then continued 
by a 1-wave from U2 to U3, followed by a 2-wave from U3 to Ur. Thus, the solution is given by 
the formula 

Ri{Ul,Ui) ® W^{Ui,U2) ® Wi{U2, C/3) ® T4^2(C/3, Ur). (3.9) 

See Figure [Tj The construction makes sense if Ai(C/2, U^) > 0, which means U3 has to be above 
U^ on yVi(C/2)- This construction can also be extended if {Ur) lies entirely above Wi(C/2). 
In this case let / and J be the intersection points of Wi{U2) and yV^(f/fl') with the axis {h — 0}, 
respectively: 

{/} ^ WiiU2) n{h^ 0}, {J} = wHUr) n{h^ o}. 

Then, the solution can be seen as containing a dry part Wo{I,J) between / and J. Thus, the 
solution in this case is 

Ri{Ul,Ui) © WsiUi, U2) © VFi([/2, /) © Woil, J) © i?2(J, Ur). (3.10) 



If ol < Ur a solution of another type can be constructed as follows. To each U G C+, a 
stationary contact to U° G G2 downing a — or to a — gl is possible, since gr > a^. The set of 
all these states U° form a curve denoted by C^. Let 

{c/i} = Wi(c/L)nc^. 
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Figure 7: Construction (Bl), > ajj; a solution of the form l|3.9 



Then, the solution begms by a 1-wave Wi{Ul, Ui), foUowed by a stationary jump W3{Ui, U2 ~ U°) 
to U2 e C+. Let {J/s} = >Vi(J72) n (?7_r). The solution is then continued by a 1-rarefaction 
wave from U2 to C/3, followed by a 2- wave from C/a to Ur. Thus, the solution is given by the 
formula 

Wi{Ul, Ui) © W^3(C/i, C/2) ® i?i((72, t/g) ® M/2(C/3, f/fl). (3.11) 

See Figure [sj The construction makes sense if \i{U^) > 0, or U3 G Gi. This construction can also 
be extended if (C/i?) lies entirely above yVi(C/2). In this case let / and J be the intersection 
points of Wi{U2) and W:f (J7i?) with the axis {h — 0}, respectively: 

{/} = WiiU2) n {h - 0}, {J} = wHUr) n{h^ o}. 

Then, the solution can be seen as containing a dry part Wo{I,J) between / and J. Thus, the 
solution in this case is 



WiiUL, Ui) © W3{Ui,U2) © RliU2,I) © Wo{I, J) © R2{J, Ur). 



(3.12) 



The wave structure of the solutions (3.9) and (3.11) are the same, but the state at which the 
solution reaches the strictly hyperbolic boundary C using a different wave. However, one may 
argue that in both cases the solution uses a stationary contact to reach C+ from either side of 
C+. Moreover, all the states in the solution {Ul, C/i, C/2, fJs, Ur) can be in an arbitrarily small ball 
center on C+ . 

Construction (B2). This case holds only when > or. Again, there is an interesting phenomenon 
as wave speeds associate with different characteristic fields coincide and all equal zero. The solution 
therefore contain three waves with the same zero speed. 

The solution begins with a 1-rarefaction wave until it reached C_|_ at Ui. At [/i, the solution may 
jump to Gi using a "half-way" stationary wave to a state M = M{a) — U°{a) from the bottom 
height ol to any a G [aR,aL]- Then, the solution can continue by a 1-shock with zero speed from 
M to TV = N{a) = N*{a) G , followed by a stationary wave from to P = P{a) = N°{a) 
with a shift in a-component from a to qr. The set of these states P{a) form a curve pattern 
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Figure 8: Construction (Bl), < a^: a solution of the form jsTTTJ. 

L. So, whenever ^ (?7h) n £ = {P = P(a)}, there is a Riemann solution containing three 
zero-speed waves of the form 

i?i(C/L, C/i) © M^3((^i, ® 5i(M(a),iV(a)) ® M^3(A^(a), P(a)) ® M^2(-P(a), t^i?). (3.13) 

See Figure |9] 

Observe that this solution coincides with the one in Construction (Bl) if the first stationary 
wave from Vi to M — U2 shifts a-component from directly to an.. The other limit case where 
the first stationary wave to Gi is not used gives a connection to the following possibility. 
Construction (B3). In this case the Riemann data can be altogether in a arbitrarily small ball in 
G2- Assume first that ol > an. Let 

Ui ^ Wi{Ul) n C+, and G G2 + resulted by W:i{Ui,U1), 

K = Wl{UL)r^C^, and A"" e G2 - resulted by W3{K,K°). (3.14) 



From any state U G Wi(C/l), where Xi{Ul) < {Ul is below Ui or coincides with Ui), we use a 
stationary jump to a state U°, shifting the bottom height from down to an. The set of these 
states U° form a "composite" curve Wi (C/l) as defined by ( |3.4[ ). The curve Wi (J/l) is thus a 
path between and K°. Whenever 7^ Wf ([7^) n (C/l) = a Riemann solution can 

be determined by 

w,iUL,UM) © WsiUM, K,) e W2{uIj, Ur), (3.15) 

where Um & >Vi(f7L), provided C/^ € G2 or A2(C/Xf > t^i?) > 0- See Figure[l0) 

Second, consider the case < a^. Let C" as in the case for the solution of the type (3.15). 



To each U G C±, a stationary contact to U° G G2 downing back a = an to a = is possible, 
since an > ql- The set of all these states U° form two curves denoted by CJ. Let 

{C/i}- Wi(C/L)nC;, and € C+ resulted by W3{U?,Ui), , . 

KeWiiUL), and K° e C- resulted by W3{K°,K) ^ ' 

decreasing an to ol- From any state U € yVi{UL), Ai(J7) < Xi{Ui), there is a stationary jump to 
a state J7°, shifting the bottom height from to a_R. The set of these states U° form a composite 
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Figure 9: Construction (B2), > aji: a solution of the form JsTTs} 
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Figure 11: Construction (B3), aj^ < ajj; a solution of the form |3.17[ ). 



curve also denoted by (C/l). Whenever ^ (f/^) n W5'(C/l) = {J7|:^}, a Riemann solution 
can be determined by 



W^(Ul, Um) ® W^{Um. Uli) ® W2{Uli, Ur), 
where Uai G Wi(C/l), provided Ur G G2 or hiU^j, Ur) > 0. 



(3.17) 



Remark 5. In both cases > aR and < qr, the condition for yV'^(C/fl) n (C/i) 7^ 
that Uf lies above Wj^Ci/ij) and K" lies below >V2^(;7ii). See Figure 



n 

Let us now discuss the existence and uniqueness. Assume first that < aR. In this case, only 
Constructions (Bl) and (B3) are available. The limit case of ( |3.9[ ) of (Bl) when = U2 coincides 
with the limit case of (3.15) of (B3). Thus, the union yVi(C/2) U 'Wi{Ul) form a continuous 
decreasing curve (the curve can be considered as the graph of u being a decreasing function of h) 
and that >Vi(f72) and (J7l) meets only at one point U2. Since yV^(?7_R) is an increasing curve, 
there always a unique intersection point of (Ur) and the union Wi(t/2) U W^j U L) if K° lies 
below or on the curve {Ur). This implies that the Riemann problem for (1.1)-(1.2) always 
admits a unique solution if K° lies below or on the curve 1^2 (Ur). 

Next, assume that ql > aR. Let U° € G2 denote the state resulted from a stationary wave 
from Ui G C+. Observe that both Ui and U* belong toWaiUi). Whenever U° lies above U* on 
>V3(C/i), there are three distinct solutions. Otherwise, there is at most one solution. See Figure 

\m 

Theorem 3.2 (Riemann problem for the shallow water equations). Given a left-hand state Ul G 
G2. 



(a) Existence. The Riemann problem (1.1)-[1.3) admits a solution if K° lies below or on the 
curve W^{Ur), and that ifW^UR) intersects with ([/l) at some point U^ G G2 , then 
A2 ([/£/, f/fl) >0. 



(b) Regime of uniqueness. The Riemann problem (1.1)-(1.3) has at most one solution if 
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either < ur; 

or > CLfi, hi > h1^ , where U2 is defined in (3.9); 



- or > aji, hi < hl^, and the states U° and Uf are located on the same side with respect 
to the curve (J/^). 



(c) Multiple solutions. If ol > clr., hi < hf, and U° lies abov e the c urve WHUr) and U. 



# 



lies below the curve W;^([/fl), then the Riemann problem {1.1)-(1.3) has three solutions. 
Example. W e provide some numerical experiments computing h1,hf to illustrate the cases of 



Theorem 



3.2 



(see Tables Bl, B2, and B3). All experiments give the same result: hi > /i* 



Table Bl 


States 


Ul 


u? 


u* 


Water Height h 
Velocity u 
Bottom Level a 


3 

0.5 
1.1 


1.819500899801235 
3.032474262659020 
1.000000000000000 


1.768961248574716 
3.119112786658156 
1.000000000000000 




Table B2 


States 


Ul 


u? 


u* 


Water Height h 
Velocity u 
Bottom Level a 


3 

0.1 
1.1 


1.707571536932233 
2.901359698616083 
1.000000000000000 


1.656818524474798 
2.990236508448978 
1.000000000000000 




Table B3 


States 


Ul 


u? 


u* 


Water Height h 
Velocity u 
Bottom Level a 


3 
1 
2 


3.187878980786353 
1.969891931767155 
1.000000000000000 


2.574902018055705 
2.438841182952260 
1.000000000000000 
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Remark 6. We conjecture that hi > h* . If this conjecture holds, then Theorem 3.2 miphes that 
the Riemann problem always has at most one solution for Ul G G2. 

3.3. Remark on the continuous dependence of solutions 

As seen in the previous subsections, the construction of Riemann solutions is based on a 



given left-hand state Ul. The Riemann problem for (1.1)- (1.2) admits up to three solutions for 
data in certain regions, which implies that the initial-value problem for (1.1)-(1.2) is ill-posed. 
However, connectivity between the types of Riemann solutions helps to determine the continuous 
dependence of the set of solutions on Riemann data. Observe that for each construction (A) and 
(B), the general structure of solutions changes continuously when changes and one evolves from 
one case to another. For example. Construction (Al) changes continuously to (A2), and (A2) itself 
changes continuously to (A3). Similar remarks hold for the cases (Bl), (B2), and (B3), as observed 
earlier about the continuity of the wave patterns. Thus, the set of solutions "globally" depends 
continuously on the right-hand side Ur for each case J7i, € Gi U C+ and Ul G G2. In order to 
show that the set of solutions depends continuously on Riemann data, we need only to check that 
when Ul moves from Gi to G2, the change in the structure of solutions is continuous as well. But 
this fact also holds true, since when U l tends to C+ on each side, the solutions in Constructions 
(Al) and (Bl) approach each other, and the solution of Constructions (A3) and (B3) approach 
each other as long as the solutions make sense. If ol > clji, these solutions eventually coincide on 



4. A Godunov-type algorithm 

4-.1. A well-balanced quasi- conservative scheme 

Given a uniform time step At, and a spatial mesh size Aa;, setting Xi — iAx,i £ Z, and 
t" — nAt,n £ N, we denote C/" to be an approximation of the exact value U{xi,t"'). Set 



The system (1.1)-(1.2) can be written in the compact form 

dtU + d^F{U) = S{U)d^a, t>0,x€R. (4.1) 
Let us be given the initial condition 

U{x,0) = Uo{x), xeR, (4.2) 

Then, the discrete initial values Uf are given by 

1 /■^i+1/2 

U° = — Uo{x)dx. (4.3) 



Ax 



-1/2 



Suppose [/" is known and [/" is constant on each interval {xi^i/2TXi^i/2) vi i e Z. On each cell 
i) we determine the exact solution to the Riemann problem for 

dtU{x,t) + d^F{U{x,t)) = S{U)d^a, on i? x (i", T+i], (4.4) 

subject to the initial condition 

C/(x,n-( S^' /^T'^''"' (4-5) 

Denote this solution by U{x,t; U"_i, U"). Use these solutions of the local Riemann problems we 
define the function V by 

V(.f\■^^ U{x,t;U:U,U:^), x,^y2<x<x,,t^ <t<t^+\ 
''y^^^i- \ U{x,t;Ur,Ur+i), x,<x<x,+i/2,t^ <t<t"+\ 
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As for the initial values, we have to ensure that the approxiniation f/"'*'^ at time t"'~^^ is constant 
on (xj_i/2: •'2^4+1/2) for a-U « G Z. Therefore, we define the new value U"^^ at the time t = by 



1 
Ax 



i+1/2 



(4.6) 



This means U^~^^ is the mean value of F on (xi_i/2, a;i_(-i/2) and thus contains parts oiU{x,t; U^_i, U^) 
and U {x, t; U^, C^i+i)- To ensure that the solutions of two consecutive local Riemann problems do 
not coincide, we assume that the following CFL (Courant, Friedrichs, Lewy) condition holds: 

At , 1 

- — max A, < -, 

where denote the eigenvalues of DF{U). 

Suppose now V is an exact solution on {xi_i/2TXi^i/2)- Since the a-component is constant in 
(2^4- 1/2 J 2:^4+1/2)1 the right-hand side of (1.1 1 vanishes for V. Thus, the standard Godunov scheme 
is in quasi- conservative form: 



up 



At 
Ax 



{F{U{x,+,/2-, t-+'; up, UP,J) - F{U{x,_,/2+, t"+'; Up_„UP))). (4.7) 



One might think that in the scheme (4.7) the source term is incorporated into the local Riemann 
problem. 



The Godunov scheme (4.7) is capable of capturing exactly equilibria. Therefore (4.7) is a well 



balanced scheme. In fact, let us be given the initial data U^ to be equilibrium states of a stationary 
wave. Then, on each cell Xj_i/2 < x < Xi^if2,t" < t < t"+^ the exact Riemann solution is 
constant. Thus, C/(a;,+i/2-, t"; C/f , ?7^i) = ?7(.Ti_i/2+, t"; C/f.i, C/f ) and so UP+^ = Up = Uf for 
all i e Z and n > 0. When there are multiple Riemann solutions, any of them can be selected and 
we still obtain a deterministic scheme, according to Theorem |3.1[ 

4-2. Numerical Riemann solver 

Given any Riemann data {Ul,Ujj), denote by U{x,t;UL,Ufi) the Riemann solution corre- 
sponding to the Riemann data {Ul, Ur). To build the Godunov scheme (4.7) we will specify the 
values J7(0±, At; Ul, Ur) for an arbitrary and fixed number At > 0. 



Riemann solver (Al). We present a computing strategy for Riemann solutions (3.1) as follows 



(i) The state U^ — {h°^,u°^,aR): — h{aB) — /ii(afl) , where /ii(afl) is the smaller root of 
the nonlinear equation (2.13), described by Lemma 2.1 and can be computed using Lemma 
= ULhL/hl. 



2.4 



(ii) The state Um 

WHUr), see ([2I0I) 



= (h M,UM,o-R) is the intersection point of the wave curves yVi(C/£) and 
Equating the it-component for these two curves leads to a strictly 
Thus, the /i-component of the intersection 



increasing and strictly convex function in h 
point /i2 can be computed using the Newton's method. 

The Riemann solver (Al) relying on Construction (Al) yields 

UiO~,At;UL,UR)^UL, 
U{Q+,At;UL,UR) = Ul. 



This implies that the Godunov scheme \A.7\ using the Riemann solver 1 becomes 

At 



UP+^ = UP 



A: 



-{F{UP)-F{{UP^,r)), 



(4.8) 



(4.9) 



where U° defined as in (4.8) 
Riemann solver (A2). 



The states of the Riemann solution (3.3) can be found as follows 



21 



(1) The state Um ~ ihM,UM,aj^) is determined by 

{UM} = m{UL)nwiiUR). 

(2) The states Ui = {hi,ui,ai),U2 — (/i2,'U2,ai) are determined by using the corresponding 
"half-way" shifting in a component from the stationary contact from Ul to U-\ an d the 
stationary contact from U2 to Um, and using the fact that U2 = Uf, (see Lemma 2.1): 



ai = aL 



Ul 



2ff 



hi — hi — Um 



2s 



+ h 



M 



_ -hi+y/hj+Shiuf/g 
tl2 - o , 



(4.10) 



"2 - h2 ■ 



solver (A2) relying Construction (A2) gives 



It is not difficult to check that the system (|4.10|) can yield a scalar equation for hi. The Riemann 

(4.11) 



UiO-,At;UL,UR)^UL, 
C/(0+, At; Ul, Ur) = Um = Um{Ul, Ur). 



This implies that the Godunov scheme ( |4.7[ ) using the Riemann solver 2 becomes 

At 



U, 



n+l 



ur - f^imn - f{Um{u-_i, um, 



(4.12) 



where UM{Ui_i,Ui)) is defined as in (4.11), i.e. 



{UMiu^-i, um = m{uiu) n wf (c/r). 

Since Um plays a key role in this Riemann solver, we sketch a computing algorithm for Um as 
follows. First we observe that if Ur lies below the curve W3{Ul) in the (/i, u)-plane, then Um is 
the intersection point of WziUL) and S2 {Ur). Otherwise, Um is the intersection point of yVziUL) 
and TZ^{Ur). Thus, we find: 

(i) (Arrival by a 2-shock) If hRUR — hLUL < then hj^i is the root of the equation 

hLUL 



0. 



(ii) (Arrival by a 2-rarefaction wave) Otherwise, /im is the root of the equation 

hLUL 



G2{h) 



{uR + 2y^{Vh- ^/hR)) =0. 



(4.13) 



(4.14) 



It is easy to see that both functions Gi,G2 defined by (4.13) and (4.14) are strictly convex. 
Moreover, we have 



G'lih) 

G'2{h) 



h2 



2 \ h 



^ < 



for all /i > 0. Thus, the Newton method can be applied for both equations (4.13) and (4.14) with 
any starting point. 
Riemann solver (A3). 

Let us consider Construction (A3) and let 



A^U 



L ' 



{B}^Wi{UL)r\W^{UR). 
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It is easy to see that J7m = (/im, "m, ol) lies on yVi([/L) between yl and _B. We propose a procedure 
similar to the Bisection method to compute the states of the elementary waves of the Riemann 
solution (3.6) as follows. We use the equation of ([/_r) : $2(C/; Ur) = 0, defined by (2.10), as 
a test condition: for U above (C/_r), ^2{U; Ur) > and for U below it, ^2{U; Ur) < 0. Using 
a stationary jump from any state U on the wave pattern of WiiUL) between A and _B to a state 
U° shifting a from to qr. Then, we have 

^2iA;UR)-<i>2{B;UR) <0. 

Algorithm 1: 

Step 1: An estimate for Hm is given by 

hA + hs 
"Af — , 



Um = {hM,UM,aL) G Wi(C/i), SO um is computed using the equation (2.6 1 with Uq = U^. 
Step 2: 

(a) If $2(v4; Ur) ■ $2(t^j\/; Ur) < 0, then set B = Um and return to Step 1; 

(b) If $2(^; Ur) ■ ^2{Um] Ur) > 0, then set A Um and return to Step 1; 

(c) If ^2{A; Ur) ■ ^2{Um', Ur) — 0, terminate the computation. 

We can still use an alternative algorithm using the value of a-component as a convergence 
condition, as follows. 
Algorithm 2: 

Step 1: Let A = Uf and B is the intersection point of Wi(?7l) and {u = 0}. An estimate for 
Hm is given by 

hA + llB 

i^M — , 



and is estimated using the equation (2.6), so an estimate of Um is Um — {^m tUm tO-l) & 
Wi(?7l). An estimate for Um is given by 

{u°M} = miUM)nwiiUR). 



Determine the change in a-component for the stationary wave between Um and J7^^ (see (2.12)) 

a — fli H z h iiM — n-M- 

Step 2: 

(a) li a — Ur < 0, then set Ha — fiM and return to Step 1; 

(b) If a — a/f > 0, then set Hr ^ fiM and return to Step 1; 

(c) If a — fl/j — 0, stop the computation. 

The Riemann solver (A3) relying on Construction (A3) yields 

U{0-,At;UL,UR) = Um = Um {Ul.Ur), 

C/(0+, At; Ul,Ur) = UIj = {Um{Ul, Ur))". ^^■'^> 



This implies that the Godunov scheme (4.7 1 using the Riemann solver 3 becomes 

Ur+' = C/« - ^{F{UM{Ur, U-^,)) F{UUU:_,, C/f))), (4.16) 



23 



where UniUr, U^+i) and Uli{Up_i,Up)) are defined as in ( [ils] ). 

Riemann solver (Bl). The Riemann solver (Bl) relying on Construction (Bl) gives 



U{G-,M-Ul,Ur) = Ui 



' UL I 2 

>3yg ^ 3 



ghL,aL t/L,+ G C+, 



If ai > a/j, then 



(4.17) 



U2 := Ul,+o e Gi, 



where Ul,+o G Gi is the state resulted by a stationary contact from Ul.+ G C+. This implies that 

L solver (Bl) becomes 

{F{U-+)~F(U-_,.,)). (4.18) 



the Godunov scheme (4.7) using the Riemann solver (Bl) becomes 

At 



n+l 



ur 



If fli < an, then U2 G C+. The computing of Ui and U2 can be done similarly as in the Riemann 
solver (A3). 

Riemann solver (B2). The Riemann solver (B2) relying Construction (B2) gives 



(4.19) 



C/(0-, M- Ul, Ur) = Ui = + yiiL) , + jVghE, flL j := Ul,+ € C+, 

U{0+, At; Ul, Ur) =P = P{Ul, Ur) e W^Ur) D W3{Ul^+). 



This implies that the Godunov scheme ( |4.7[ ) using the Riemann solver 2 becomes 

At 



jjn+l ^ ^r. 



Ax 



{F{U-^) ~ F{P{u:_„U:))), 



(4.20) 



where 



p(c/r„i, c/D) - (c/D n W3(f/;^i,+). 

Riemann solver (B3). The Riemann solver (B3) relying on Construction (B3) yields 

U{0-, At; Ul,Ur) ^Um = Um{Ul. Ur), 
U{0+, At; Ul,Ur) = UIj = {Um{Ul,Ur))°. 



This implies that the Godunov scheme (4.7 1 using the Riemann solver 3 becomes 

At 



jjn+l ^ 



Ax 



{f{Um{u:, c/i;i)) - F{UUU-_„ C/f ))), 



(4.21) 



(4.22) 



where Um{UJ^ ,Ul\^^) and Ulj{U"_i,U-'-)) are defined as in (4.21). It is easy to see that the 
determinations of the states U {Q— , At;U l ,U r) and C/(0+, At; [/l, CZ/j) by Solvers (A3) and (B3) 
are the same. 

The computing strategies for Solvers (Bl), (B2), and (B3) are similar to those of Solvers (Al), 
(A2), and (A3). 

4-. 3. The Godunov algorithm 

It is natural to ask which Riemann solvers should be taken in the Godunov scheme. As seen 
earlier, in the regions where there are possibly multiple solutions one can select any Riemann 
solution. Three extreme cases can be distinguished by preferring one particular Riemann solver 
whenever it is available. For example, we can decide to always select solutions with stationary 
contact wave in the same region as the left-hand state, that is, the solvers (Al) and (B3). This 
selection leads us to a deterministic algorithm for designing a corresponding Godunov scheme, as 
now described. 
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Building Godunov Scheme Algorithm preferring solvers (Al) and (B3). Let Ul — and — 

Tjn 

If Ai([/l) < 

If $2(f/i*;C/fl) <0 
Use Solver (Al) 

elseif $2(C/f°;C/i?,) < 

Use Solver (A2) 
else 

Use Solver (A3) 
end 
else 

If $2(C/i°;(7j^) >0 

Use Solver (B3) 
elseif $2(C/f ;t/i?) > 

Use Solver (B2) 
else 

Use Solver (Bl) 
end 
end 



5. Numerical experiments (I). The non-resonant regime 

We now numerically investigate our Riemann solver and Godunov method and present several 
numerical tests. For each test we consider the errors between the exact Riemann solution and 



the approximate solution by the Godunov scheme (4.7 1 for x e [—1,1] with different mesh sizes 
corresponding to 500, 1000, 2000 points. In this section as well as in the next section, we plot the 
solution at the time t = 0.1, and use the stability condition 

CFL = 0.75. 

The algorithm for selecting the Riemann solvers is the one described at the end of the last section, 
unless indicated otherwise. 

5.1. Test 1 

This test indicates that the Godunov method is capable of maintaining equilibrium states. Let 

^ ' \ Ur = {hR,UR,aR), X > 0, ^ ' 

where Ul = (1,5, 1) and Ur = (1.223655890827479,4.086116070277590, 1.2). It is not difficuh to 



check that the Riemann problem with initial data (5.1 ) admits a stationary contact between these 
equilibrium states: 

U{x,t) ^Uoix), xeR,t>0. (5.2) 



Figure 13 shows that the stationary contact is well captured by Godunov method using our exact 



Riemann solver for x Cz [—1,1] with 500 mesh points and at time t — 0.1. 
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Godunov scheme: water height Time =0.10012 





500 mesh points . 

Exact solution 













-1 -O.B -0.6 -0.4 -0.2 0.2 0.4 0.6 0.6 1 



Velocity. Error =0.0030874 



5.5 





500 mesh points 

Exact solution U 
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Figure 13: Test 1: A stationary contact wave is captured exactly by Godunov method using our exact Riemann 
solver. 



5.2. Test 2 

We now approximate a non-stationary Riemann solution with data Ul,Uij G Gi. Precisely, 



we consider the Riemann problem (|1.1|)-(|1.3P with data 



Ul = {hL,UL,aL) = (0.3, 2, 1.1) e d, a; < 0, 
Ur = {hR, UR, ur) = (0.4, 2.2, 1) e Gi, x > 0. 



(5.3) 



The Riemann problem (1.1 1 (1.2 1 with the initial data (5.3) admits the solution described by 
Construction 1, where 

Ul = (0.21815897, 2.750288, 1), U2 = (0.35252714, 1.9572394, 1). 

The errors for Test 2 are reported in the following table 



N 


\\U^-U\\l^ 


500 
1000 
2000 


0.012644 
0.0087928 
0.0063773 



and Figures [14] and [15] and the table above show that approximate solutions are closer to the exact 
solution when the mesh size gets smaller. 

5.3. Test 3 

In this test, we will approximate a non-stationary Riemann solution with Riemann data 



Uli Ur € G2. Precisely, we consider the Riemann problem (1.1)-(1.3| with data 

rrr^f f/L = (/iL, ML, fli) = (1,3,1.2) e G2, 2^ < 0, 

I UR = {hR,UR,aR)^ (2,0.5,1) eG2, x > 0. 
This problem admits the solution described by Construction (A3), where 

Um = (1.8452179, 0.67672469, 1.2), U^j = (2.0496463, 0.60922927, 1). 



(5.4) 
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Water height 



- Exact solution 
-500 mesh points 
■1000 mesh points 
-2000 mesh points 
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Figure 14: Test 2. Water height with different mesh sizes. 



Water velocity 
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Figure 15: Test 2. Water velocity with different mesh sizes. 
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Water height 



- Exact solution 
-500 mesh points 
■1000 mesh points 
-2000 mesh points 
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Figure 16: Test 3. Water height with different mesh sizes. 

Precisely, the solution starts with a 1-shock from Ul to C/m, followed by a stationary contact from 
Um to , and then arrived at Ur from W^j- by a 2-shock wave. The errors for Test 3 are reported 
in the following table: 





\\U^-U\\l^ 


500 
1000 
2000 


0.01813 
0.0076434 
0.0035277 



(5.5) 



Figures [16] and [17] and the table above show that approximate solutions are closer to the exact 
solution when the mesh size gets smaller. All of our tests presented so far exhibit a convergence 
of approximate solutions to the exact solution. 



6. Numerical experiments (II). The resonance regime 

In the following, we will consider the cases where the Riemann data on the different sides with 
respect to C+ : 

(i) Ul e Gi and Ur e G2; 

(ii) Ul e G2 and Ur e d. 

The solution is evaluated for x E [—1,1] with 500 points and at time t = 0.1. We take also 

CFL = 0.75. 

6.1. Test 4 

In this test, > ur, and there is a unique solution. We consider the Riemann problem 



([LlJ)-([L3j) with data 
Uo{x) 



C/l = (/iL,UL,aL) = (1,3,1.1) e Gi, x<0, 
Ur = {hR, UR, ur) = (1.2, 0.1, 1) e G2, x> 0. 



(6.1) 



We have 



= 1.042865405801653 < hf" = 1.213385283426733. 
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Water velocity 



2.5 



1.S 



D.5 



- Exact solution 
-500 mesh points 
■1000 mesh points 
-2000 mesh points 
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Figure 17: Test 3. Water velocity with different mesh sizes. 



Thus, the problem (1.1), (1.2), (6.1) admits a unique solution of the form (3.1), according to 



Theorem (3.1), where 

Um = (1.5521168, 1.4328264, 1.1), J/^^ = (1.665941, 1.3349296, 1). 



Figures [T8][T9| show that the Godunov scheme gives good approximate solutions to the exact 
solution in this resonance case. 

6.2. Test 5 

In this test, aj^ < a^i, and there is a unique solution. We consider the Riemann problem 



([rT])-([r3| with data 
Uoix) 

and 



C/l = (/iL,WL,aL) = (0.2,4,1) e Gi, a: < 0, 
Ur = {hR,UR,aB,) = (0.5,1.5,1.1) G Ga, x> 0, 



U°* = (0.677264819960833,1.181221844722815), 
^2{UI*;Ur) = -1.050411375011095 < 0, 
U*° = (0.581828763814630, 1.374974992221044), 
$2(C/f °; Ur) = -0.474326705580410 < 0. 



(6.2) 



This Riemann problem admits a unique solution of the form (3.1), according to Theorem 3.1 
where 

[/£ = (0.21591647, 3.7051366, 1.1), Um = (0.56185289, 1.7661913, 1.1). 



Figures 20 to [21] show that the Godunov scheme gives good approximate solutions to the exact 
solution in this resonance case. 

6.3. Test 6 

Next, we provide a test for the case there are multiple solutions. So we consider the Riemann 



problem (1.1)- (1.3) with data 



Uo{x) 



t/L = (/iL,UL,aL) = (0.2,5,1) e Gi, a;<0, 

Ur = (hR, UR, ur) = (0.75904946, 1.3410741, 1.2) e G2, a; > 0. 



(6.3) 



29 



Water height 



- Exact solution 
-Approximate 
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Figure 18: Test 4 (Resonant case). Water height of the solution (jSJ 



Water velocity 
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- Exact solution 
-Approximate 
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Figure 19: Test 4 (Resonant case). Water velocity of the solution (|3.6| 
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Figure 20: Test 5 (Resonant case). Water height of the solution ( |3.1| . 



Water velooity 
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Figure 21: Test 5 (Resonant case). Water velocity of the solution ( |3.1| . 
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Water height 
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Figure 22: Test 6 (Resonant case - multiple solutions). Water height of the solution l |3.1| l - preferred Solver (Al). 
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Figure 23: Test 6 (Resonant case - multiple solutions). Water velocity of the solution | |3.1| - preferred Solver (Al). 
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Water height 



- Exact solution 

- Solution by Solver 2 
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Figure 24: Test 6 (Resonant case - multiple solutions). Water height of the solution ||3.3| preferred Solver (A2) 



The Riemann problem (1.1)-(1.3) with mitial data (6.3) admits three distinct sohitions: one 
solution of the form (3.1) with 

Ul = (0.21984063, 4.5487497, 1.2), Um ^ (0.7964266, 1.4737915, 1.2), (6.4) 



one solution of the form (3.3) with 

Um = (0.75904946, 1.3174372, 1.2), 



(6.5) 



which can be seen to be a stationary solution, and one solution of the form ( |3.6[ ) with 

Um ^ (0.95328169, 0.89892673, 1), U^j = (0.72279573, 1.1855776, 1.2). (6.6) 

We could have three extreme choices of a Riemann solver for the Godunov method in this case. 
This can be seen easily by saying that we prefer a particular solver whenever it is available. From 



Figures 22 and 23 it follows that if Solver (Al) is preferred whenever it is available, then the 
approximate solution approaches the exact solution. The same observation for Solver (A2); see 



Figures 24 and 25 However, it is not the case for Solver (A3): approximate solutions do not 



converge to the exact Riemann solution by (3.6); see Figures 26 and 27 



6.4. Test 7 

Finally, we consider the Riemann problem with data 



Uo{x) 



We have 



C/l = (/JL, ML, fli) = (1,2,1.1) e G2, a;<0, 
Ur = {Hr, ur, ur) = (0.8, 4, 1) e Gi, x> 0. 

0.998204556070240 <hl^ 1.050890579855180, 



(6.7) 



where hf,h° are defined as in Theorem (3.2). Thus, the problem (1.1 1, (1.2 1, (6.1) admits a 
unique solution of the form (3.10), according to Theorem (3.2), where 



Ul = (0.77374106,2.7536634, 1.1), C/2 = (0.58589019,3.636556,1), 
U3 = (0.64142927,3.4143821, 1). 



(6.8) 
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Figure 25: Test 6 (Resonant case - multiple solutions). Water velocity of the solution (|3.3| - preferred Solver (A2) 
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Figure 26: Test (Resonant case - multiple solutions). Water height of the solution \3.G\ - preferred Solver (A3) 
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Water Velocity 




Figure 27: Test 6 (Resonant case - multiple solutions). Water height of the solution \3.6\ - preferred Solver (A3) 



Water height 




Figure 28: Test 7 (Resonant case). Water height of the solution l |3.10| . The Godunov scheme does not approach 
the exact solution in this resonance case. 
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Water velocity 




Exact solution 

ApproM. solution 
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Figure 29: Test 7 (Resonant case. Water velocity for | |3.10[ l. The Godunov scheme docs not approach the exact 
solution in this resonance case. 



Precisely, the solution begins with a 1-rarefaction wave from Ul to C/i, followed by a stationary 
contact from Ui to J72, then continued by a 1-shock wave from U2 to 6^3, and finally attains Un 
by a 2-rarefaction from t/3. Figures [28p9| show that the Godunov scheme generates approximate 
solutions which do not approach the exact solution in this resonance case (here shown at time 
t = 0.03). 

7. Concluding remarks 

In this paper, we have provided for the first time a complete characterization of all solutions to 
the Riemann problem associated with the shallow water equations. First, we determined domains 
in the phase space in which precisely one solution exists and domains in which several solutions (up 
to three) are available. Second, we provided a computing strategy which allowed us to numerically 
determine all possible Riemann solutions. Third, we defined a well-balanced and quasi-conservative 
Godunov scheme, which was carefully tested in several regimes of interest. 

The following main conclusions were established in this paper: 

• The proposed scheme captures exactly the equilibrium states. 

• It converges to the uniquely defined solution to the Riemann problem in strictly hyper- 
bolic domains of the phase space, and this property validates our numerical strategy. We 
emphasize that the existing literature restricts attention to this strictly hyperbolic regime, 
only. 

• Next, in order to further test our new algorithm, we considered Riemann data belonging to 
both sides of the resonance curve; both convergence to the selected solution and as well as 
convergence to another solution were observed. This is not surprising since multiple solutions 
are available in such a regime. 

• Finally, the most challenging test was performed by taking states precisely on the resonance 
curve, and we observed that the scheme gave quite unsatisfactory results with no convergence 
whatsoever observed. This latter behavior should not be interpreted as a drawback of the 
numerical method itself, but rather indicates a limitation of the physical model which, in 
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itself, does not properly describe the fluid flow so that further physics is required in this 
regime. It is conceivable that a more satisfactory model could be obtained by analyzing the 
(small-scale) effects of higher-order terms and, for instance, introducing a suitable notion of 
kinetic relation, similarly to what is done for undercompressive shock waves [251 
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